The thermomagnetic instability in superconducting films with adjacent metal layer 
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Dendritic flux avalanches is a frequently encountered consequence of the thermomagnetic instabil- 
ity in type-II superconducting films. The avalanches, potentially harmful for superconductor-based 
devices, can be suppressed by an adjacent normal metal layer, even when the two layers are not 
in thermal contact. The suppression of the avalanches in this case is due to so-called magnetic 
braking, caused by eddy currents generated in the metal layer by propagating magnetic flux. We 
develop a theory of magnetic braking by analyzing coupled electrodynamics and heat flow in a 
superconductor- normal metal bilayer. The equations are solved by linearization and by numerical 
simulation of the avalanche dynamics. We find that in an uncoated superconductor, even a uniform 
thermomagnetic instability can develop into a dendritic flux avalanche. The mechanism is that a 
small no n- uniformity caused by the electromagnetic non-locality induces a flux-flow hot spot at a 
random position. The hot spot quickly develops into a finger, which at high speeds penetrates into 
the superconductor, forming a branching structure. Magnetic braking slows the avalanches, and if 
the normal metal conductivity is sufficiently high, it can suppress the formation of the dendritic 
structure. During avalanches, the braking by the normal metal layer prevents the temperature 
from exceeding the transition temperature of the superconductor. Analytical criteria for the insta- 
bility threshold are developed using the linear stability analysis. The criteria are found to match 
quantitatively the instability onsets obtained in simulation. 

PACS numbers: 74.25.Ha, 68.60.Dv, 74.78.-w 



I. INTRODUCTION 

The concept of the critical state introduced by Bean^ 
is widely used to describe various physical properties in 
the vortex phase of type-II superconductors, see, e.g., 
Refs.^and J, and references therein. According to Bean, 
the driving force of the currents is balanced by the pin- 
ning force from material inhomogeneities, with a strength 
characterized by the critical current density, jc- The dis- 
sipation becomes vanishing if the current density j is less 
than jc. 

However, the critical state can be unstable with re- 
spect to fluctuations, e.g., in the temperature. Since 
jc decreases with increasing temperature, a fluctuation 
causing an increase of the temperature will facilitate 
further penetration of magnetic flux into the sample, 
and consequently more dissipation. This positive feed- 
back loop is the mechanism behind the thermomagnetic 
instability,'^ see also Refs. [6| and [7l for a review. In 
bulk samples the instability often results in large flux 
jumps, sometimes causing the entire superconductor to 
heat above its superconducting transition temperature, 
Tc^ In thin films, the instability leads to flux avalanches 
showing complex dendritic structures.^ ^^ In both cases 
the system is bistable since it, after some time, reaches 
a stable highly dissipative state characterized by linear 
response of the current to electric field.^ Criteria for 
onset of the thermomagnetic instability were first con- 
sidered for bulks under adiabatic conditions.^'Sn^ The 
theory was later extended to include also the fiow of 
heat ,^^^H^ and it was found that the instability onset can 



be accompanied by oscillations in temperature .^^^^^ For 
superconducting films in transverse geometry, the anal- 
ysis is complicated by the fact that the electrodynamics 
is nonlocal. ^^ ^^ Also, close to the critical state, the re- 
sponse of the electric field to fluctuations in the current 
density is strongly nonlinear .l^^^ Therefore, one has to 
address an essentially nonlocal and nonlinear dynamical 
problem. Nevertheless, linear stability analysis has suc- 
ceeded in providing criteria for the lower threshold fleld 
for onset of the instability and the upper threshold tem- 
perature, which may be signiflcantly lower than J'c.'^^^^ 
It has also been shown that edge defects can lower the 
threshold for onset of the instability.^^ ^^ The time evo- 
lution of dendritic flux avalanches has been investigated 
using numerical simulations, which produce patterns in 
striking resemblance with experiments.^^ ^^ 29 32 

To prevent the thermomagnetic instability to occur, 
superconducting cables are constructed by embedding 
superconducting fllaments in a normal-metal matrix.^ 
Indeed, experiments also show that a metallic layer is 
beneflcial for the stability of superconducting fllms, as it 
can completely suppress dendritic flux avalanches j ^^ * ^^ ' 
The suppression becomes more efflc ient a s the thick- 
ness of the metallic layer incr eases. ^^^^ It has also 
been reported that avalanches change direction when 
meeting the metal-covered part of a superconducting 
sample. ^^^^ The proposed mechanism behind the ob- 
served suppression was based on the idea that the 
normal- metal layer acts as a thermal shunt, decreas- 
ing temperature gradient s.^^ However, experiments have 
demonstrated that the avalanche activity can be reduced 
even when there is a spatial gap between the metal and 



the superconductor!^ This shows that avalanches can 
be prevented also by the eddy currents induced in the 
metal, i.e., a magnetic braking effect. Further evidences 
of the avalanche-induced non-stationary eddy currents in 
an adjacent metal layer, are the voltage pulses appear- 
ing when branches of dendritic flux avalanches take place 
under the metal layer in a partly metal-covered supercon- 
ducting film.^^ 

In the present work we will investigate the stability of 
the flux distributions in a superconducting film located 
close to a normal- metal layer. The analysis focuses on the 
process of magnetic braking, so we will assume that there 
is no thermal contact between the two layers. Magnetic 
braking is a dynamic effect, so in order to understand how 
it operates it is necessary to compare the time evolutions 
of metal-covered and pristine superconductors. Hence, a 
large portion of the present work is devoted to the study 
of pristine superconductors. 

To reach our goal we will use a combination of linear 
stability analysis and numerical simulations. The linear 
stability analysis will provide the conditions for onset of 
instability, while the simulations allow us to follow the 
full time evolution, from nucleation of the instability to 
the formation of the dendritic structures. The analysis is 
repeated with a metal layer included, but under otherwise 
identical conditions. We will thus find how the magnetic 
braking affects dendritic flux avalanches. 

In the study of the pristine superconductors, we will 
pay particular attention to the long wavelength modes of 
the Fourier space linear stability analysis. Such modes 
are particularly interesting because they become unsta- 
ble at low electric field, but they have been neglected in 
previous linear stability analysis focused on nucleation 
of finger-like patterns elongated transverse to the edge. 
A major question addressed in the present work is if a 
uniform thermomagnetic instability can develop into a 
dendritic structure. 

The paper is organized as follows. Section II describes 
the model and outlines the setup for simulations. Sec- 
tion III finds a formal solution of the equations to first 
order in the perturbations. Section IV considers the sta- 
bility of uncoated superconductors. Section V considers 
the stability of metal-coated superconductors. Section VI 
gives the summary and discussion. 



II. MODEL 

A. Basic Equations 

Let us consider a superconducting strip with an ad- 
jacent metal layer, as depicted in Fig. [l] For simplicity 
we assume that there is no thermal coupling between the 
superconductor and the normal metal, while the super- 
conductor is thermally coupled to the substrate, which 
is kept at constant temperature Tq. The thickness of 
the metal, d^, and superconductor, d^, are both much 
smaller than the strip width, 2w. Therefore, we can 
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FIG. 1. (color online) Sketch of the system: a thin super- 
conducting strip of thickness ds with a deposited metal layer 
of thickness dm- The superconductor is in thermal contact 
with the substrate, kept at constant temperature To, but not 
with the metal. Current flows in the y direction and flux has 
penetrated a distance Ix from both sides due to the applied 
magnetic field Ha. 



parametrize the problem using the sheet current J, de- 
fined as j = J(5(z), where j is the total current density 
and 6{z) is the Dirac delta function. 

The sheet current J entering the Maxwell equations 
consists of two contributions,^ 



J = Js 



(1) 



where and Jg and Jm are the sheet currents in the su- 
perconductor and metal layer, respectively. Since the two 
layers are close, the electric field, E, is approximately the 
same in the two layers, giving 



J 5 — dgCTgrj^ J 772 — (^m^m^' 



(2) 



The conductivity of the metal, <Jm, is assumed to be E- 
independent. The current- voltage relation in the super- 
conducting film is assumed to be non-O hmic with E- 
dependent conductance expressed as^^^ 



^ I {Eds/pnJc) 
Pn I 1, 



l/Us 



J < JcdiudT < Tc, 
otherwise . 



(3) 



Here T is the local temperature, Jc = djc is the sheet 
critical current, pn is the resistivity of the superconductor 
in the normal state, and n^ is the creep exponent. 

The electrodynamics is governed by the Maxwell equa- 
tions, 

VxE = -B, V-B = 0, VxH = J(5(z), (4) 

with B = /ioH and V • J = 0. The flow of heat in 
the superconductor is described by the energy balance 
equation describing the interplay between Joule heating, 
thermal conduction along the film, and heat transfer to 
the substrate. It reads as 



ct = K.V'^T 



h 



{T - To) 



- J.s • E , 



(5) 



with superconductor specific heat c, heat conductivity a>:, 
coefficient of heat transfer to substrate h. Since there is 
no thermal contact between the metal and the supercon- 
ductor there is no need to calculate the flow of heat in 
the normal metal. 



the ramp rate of the applied magnetic field, E r^ Hai^ 
This relationship is needed to bring together the linear 
stability analysis with numerical simulations and exper- 
iment. In most experiments the ramp rate is moderate, 
say Ha <C IQ-"^. 



Dimensionless form 



C. Numerical procedure 



For further analysis it is convenient to express the 
equations in a dimensionless form. We denote 

T - I_ J - A_ 7-A fr - ^ ~_^ ~_1 

Ic JcO JcO JcO ^ ^ 

7 . Pn ^ E ^ ^ dm 

t = t , rj = ^, as = (JsPm (^m = O-mPn—r ' 

jJ^odsW PnJcO ds 

Here Jco is the sheet critical current at T = 0. Henceforth 
we omit the tildes for brevity. In these units the heat 
propagation equation reads as 



f = aV^T - p(T - To) + -i^JsE, 



(6) 



where 7 = c{Tc)/c{T) is a function of temperature and 
a, /3, and 7 are constants, provided the ratios n/c and 
h/c are independent of temperature (that we assume). In 
Eq. ([6|, Of is dimensionless heat conductivity, j3 is dimen- 
sionless constant for heat transfer to the substrate, and 
7 is the Joule heating parameter. The dimensionless ma- 
terial parameters are related to the physical parameters 
as 



li^nd 
Pncw' 



P 



PnC 



!■ 



TrC ' 



(7) 



where all quantities are evaluated at Tc. 
The dimensionless Maxwell equations are 



V X E = -H, V-H = 0, VxH = JS{z) . (8) 

The material laws can be expressed in the dimensionless 
form as 



Js 






(S/Jc)i/"% J< Jc and r<l, 



otherwise , 



J, 



m — '^m 



E. 



(9) 



The above expressions are valid for arbitrary temperature 
dependencies of Jc^ Us and 7. To be specific, we will 
assume cubic temperature dependencies for /^, /i, and c, 
and linear temperature dependency for Jc^ as typical for 
low-Tc superconductors and MgB2, i.e., 

Je = l-T, n, = ni/T, 7 = T-^ (10) 

The parameters used in the calculations of this work 
are chosen to be compatible with the formation of den- 
dritic structures. For example, we let a = 10~^, /3 = 0.1, 
7 = 10, and rii = 20^ In our analysis, the electric field 
is kept as a free variable. However, one should keep 
in mind that in the critical state it is proportional to 



The simulations are performed for an infinitely long 
strip extended in the y direction, as depicted in Fig. [l] 
We analyze the full nonlinear problem by numerical 
time integration of Eq. (IgI) (the heat flow equation) and 
Eq. ([8| (the Maxwell equations) with the material re- 
lations given by Eq. ([9| and temperature dependencies 
given by Eq. (10). The set of boundary conditions and 



the calculation procedure are detailed in Ref. [29] 

In order to make the comparison with the linearized 
theory as close as possible and elucidate dynamics of the 
dendrites, the numerical analysis is conducted in two sep- 
arate steps. 

At the first step we find the background flux distri- 
bution by solving the Maxwell equations decoupled from 
the thermal effects. For that we put 7 = 0, starting from 
a zero- field-cooled state, and ramp applied magnetic field 
with a constant rate Ha until the flux has penetrated over 
a given distance Ix. 

At the second step the thermal feedback is turned on 
by putting 7 > 0. The state will then start evolving, 
and the difference from the background state is called 
the perturbation. The background state is stable if the 
perturbation saturates to some small value and unstable 
if it develops into a dendritic flux avalanche. 

It is worth noting that the formulated numerical pro- 
cedure differs from the conventional linear stability anal- 
ysis in several aspects. In particular, (i) the pertur- 
bations (as they are defined above) are not necessarily 
small; (ii) the background distributions of 5^, E, J, and 
T are essentially non-uniform; (iii) The electromagnetic 
boundary conditions are more proper; (iv) due to flux 
creep, the maximum current density in the critical state 
is slightly lower than the critical current density, i.e.. 

The numerical procedure of this work deviates 
from previous numerical simulations of dendritic flux 
avalanches in the absence of randomly distributed dis- 
order. Such disorder is important since it causes fluc- 
tuations in the background £^- values, which may trigger 
avalanches. ^^ ^^ However, for simplicity of the calcula- 
tions and the analysis, the present work considers only 
spatially uniform samples. 



III. LINEAR STABILITY ANALYSIS 

Let us assume that we start from a uniform background 
distributions of the electric field E = Ey and tempera- 
ture T, as depicted in Fig.[l] Due to the applied magnetic 



field or current, the magnetic flux front, and tiius also the 
fronts of E and T have reached a distance Ix from both 
edges. The perturbed values of E and T are specifled as 
E + JE and T-\-ST. To meet the boundary conditions we 
assume that in the Fourier space the perturbations are 
of the form 



^Ex = €x e^^ sm{kxx) sm{kyy), 

SEy = Ey e^^ cos{kxx) cos{kyy)^ 

ST = Oe ^ cos{kxx) cos{kyy)^ 



(11) 



where kx and ky are the in-plane wavevectors and A is the 
instability increment. The flux penetration depth sets 
the lower limit for allowed wave-vectors in x direction 
and we will thus identify Ix = 7T/2kx. The electrical 
current and magnetic fleld perturbations are 



6Jy 

SH, 



oAt 



SJx = ix e sin{kxx) s'm{kyy), 
cos{kxx) cos{kyy)^ 
b e^* sin{kxx) cos{kyy). 



(12) 



We will now linearize the equations in the perturba- 
tions and find a solution for the instability increment 
X = X{E,T,kx,ky). 

After linearizing the product JsE in Eq. ^ and mak- 
ing Fourier transform we express the heat propagation 
equation as 



1 11 JsE 



1 



L (13) 

where 1/T* = \d\nJc/dT\ and k = Jk'^^k'^. We have 



used that 5\'E\ = 5Ey due to the boundary condition 
E • X = 0. The temperature derivative of 7 has been 
ignored. 

The perturbations of the current components ix and 
iy are related to £ and through Eq. ^ (the material 
laws) , which gives for the Fourier amplitudes of the per- 
turbations 



J 

E' 



_J_ 



1 Js 



(14) 



Here we have introduced the nonlinearity exponent of the 
composite system, n = n{E,T), as 



d\nE 
dlnJ 






(15) 



When electric field is small most current fiows in the 
superconductor, giving n ^ ris- For high electric fields, 
the current may fiow in the normal metal and this is the 
regime in which we expect that magnetic braking can 
suppress the thermomagnetic instability. 

The value of Sx is fixed by requiring continuity of the 
current, V • (5J = 0, giving kxix 



kyiy = 0. Thus 



Eky . 

u rxix 



ky 
fVnr._ 



1 



lEJs 



JT* 



(16) 



Using the Faraday law 

Xb = kxSy -\- ky£x 



and Eq. (16), we get 



^^ = f P^ + f 



IE J, 



JT* 



(17) 



(18) 



The Biot-Savart law relates Bz with J. Treating the film 
as infinite gives the simple relation^ 



1 k 



t^yt 



(19) 



where the continuity of current has been used to elimi- 
nate ix' The exact treatment of the film boundary would 
transform the above relation to a sum over k' ^ where 
the off-diagonal elements with k' ^ k are largest for the 
longest wavelengths.^ However, for the linear stability 
of the system, the diagonal elements are by far the most 
important, and we will keep only those in the stability 
estimates. Using the expression for the current, Eq. (14), 
we get 



k 1 J 
2 kx nE 



k 1 Us — 1 Js 

Z fvx '^S J- 



(20) 



Combining Eqs. (18) and (20) we eliminate 6, and after 
some algebra we get 



Ky 



kjj^^ 

2nE' 



1 Js 



^l- 




^ A 


+ 


-X 


^ J 




'^ \ 



Combining the above equation following from electrody- 
namics with Eq. (13), describing the heat fiow, gives a 



quadratic eigenvalue equation for A, 

AX^ ^BX^C = 0, 
where 
A 

B 
C 



(21) 
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"2 V 


Us J 


n 


( A-2\ 
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ns_ 


-lljJc 









-JF, 



rji^ 



(22) 



According to the temperature dependencies (10), F = 
jT~^ when n^ ^ 1. 



The largest of the two solution of Eq. (21 ) is 



(23) 
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FIG. 2. (color online) The threshold for onset of instability in 
the T — E plane, calculated as contours of max{Re A(k)} = 0. 
The state is unstable to the left of the contours and stable to 
the right. As indicated in the figure, one of the parameters a, 
/3 , 7, /cc and rii is changed in each panel. In all panels, the 
middle curve is the same, and it is created with a = 10~^, 
13 = 0.1, 7 = 10, la: = 0.1, m = 20. 



A mode is linearly stable if Re A < and linearly unsta- 
ble if Re A > 0, and the threshold condition for instability 
is Re A = 0. If B'^ > AAC, A is real, and the threshold 
condition is C = 0. li B'^ < 4 AC, A is complex, and 
threshold condition is B = 0. The existence of an imag- 
inary part means that the state will be oscillating. The 
frequency is a; = Im A = ^/C/A. 



IV. UNCOATED SUPERCONDUCTOR 

In this section we will consider an uncoated supercon- 
ductor which can be accounted for by putting (Jm = 0. 



A. Instability threshold 

Figure [2] shows contours of the the instability thresh- 
old max{ReA(k)} = 0, with T on the horizontal and E 
on the vertical axis, obtained by numerical solution of 
Eq. ( [23| ). The instability increment value at the most 
unstable mode, max{ReA(k)}, is found by iterating over 
a finite number of ky values. The value of kx is fixed as 

The panels demonstrate how the instability threshold 
contours shift when changing a, /3, 7, ni and Ix- In all 
panels the middle curves are the same and can be used as 
a reference. It has the parameter combination a = 10~^, 
13 = 0.1, 7 = 10, Ix = 0.1, and rii = 20. 

As expected from both experiments and previous linear 
st ability a nalysis, the system is stable for low E and high 
y| i7 | 25 | 27 | rpj^g system becomes more unstable for increas- 
ing values of the Joule heating parameter 7, the creep 
exponent parameter ni, and width of the flux-penetrated 
region Ix. The system improves stability for increasing 
values of the lateral heat transport parameter a and pa- 
rameter for heat transport to the substrate /3. The graphs 
are equally sensitive to /3 at all T, but most sensitive to 
a at high T. This indicates that the instability threshold 
at low T is at longer wavelengths (smaller ky) in which 
the lateral heat diffusion is less important. The kinks 
in the curves occur at B^ = 4 AC, which is a transition 
point to the oscillatory regime. 

When cr^ = we have Jg = J, rig = n. Then the 



coefficients B and C in Eq. ( 22 ) simplify to 




(24) 



EF. 



At low T and E the state is close to what is described 
by the Bean model, therefore, we let n ::^ 1, J = Jc, and 
T = To. 

We will now derive closed expressions for the threshold 
electric field, Eth{T, kx), in three limiting cases. 

0. Uniform oscillatory instability. At low electric 
fields the onset of instability takes place in a thin layer 
along the edge, so that k^. ^ 1,^^- In this case 
B'^ < AAC and the solutions of the dispersion equation 
Re A = are oscillatory. The instability onset in this case 
corresponds to 5 = 0, or 



B 



kx akl + p 
~2 nE 



Jr 



i^x 



JrF = 0. 



(25) 



Solving this equation for E gives the threshold electric 
field 



E. 



(0) 
th 



Jr akl 



■13 



U Q-L ZjIxix 



(26) 



The physical interpretation of Eq. (26) is straightforward: 



increasing heat removal through a and P increase the 



threshold, while increasing Joule heating through F and 
nonlinearity through n decrease it. 

In the Bean model limit, n ^ oo, the threshold is in- 
dependent of a and p. This corresponds to the adiabatic 
limit: 



^Miab) ^ j^^/2. 



(27) 



This means that the sample is always stable for Ix < 

1. Finite-wavelength oscillatory instability. At higher 
temperatures and higher electric fields the instability will 
nucleate at smaller kx and the most unstable mode will 
be at finite ky. Still kx and ky are comparable in size so 
that the instability is accompanied by oscillations. This 
means that the condition for onset of instability is 5 = 
in Eq. (24). When neglecting the ky/n term, it becomes 



B = 



2nE 



ak^ 



JcP 

nE 



JcF]k^kl = 0. 



(28) 



The most unstable mode is found by the condition 
dReX/dky = 0, which gives dB/dky = 0, i.e.. 



(29) 




Elimination of ky leads to 




Kf^r 



0. 



(30) 



This equation straightforwardly gives the onset condition 
in terms of the threshold kx. However, in this work we 
focus on the threshold electric field, so the equation must 
be solved for E. Changing variables to x = {P/FnE)^^^ 
gives a cubic equation x^ + 3px + 2g = 0, with the coef- 
ficients p = {JcF)-'^/^{a/pY^^kt^^ and q = -1/2. This 
is solved with Cardano's formula x = ii+ + ii_, where 



u± 




-.1/3 



P {JcFf 



Thus, the threshold electric field for finite- wavelength os- 
cillatory instability is 



e'^ = 4^{u^ + u.)- 



(31) 



Due to the approximations used in the derivation, the ex- 
pressions is mainly of value for small kx . Series expansion 
gives 






Fn 



P J^F^ 



1/3- 



(32) 



The peculiar kx dependency is due to the k/2 Biot-Savart 
kernel and is thus a consequence of the nonlocal electro- 
dynamics. 



2. Fingering instability. When the temperature is 
sufficiently high the oscillations cease to exist and the 
instability is elongated transverse to the edge with ky ^ 
kx. It is thus called a fingering instability and the condi- 
tion for onset is C = in Eq. (24): 



C = {akl 




EFkl = 0. (33) 



The most unstable mode is at dX/dky = 0, giving 
dC/dky = 0. Hence, 

2akl = nEF - nakl - /3. (34) 

Eliminating ky and solving for E gives the threshold elec- 
tric field for the fingering instability 



e: 




(35) 



This is the same expression as found in Refs. ^24^ and |25l 
In order to compare the above threshold conditions 
with experiments, they must typically be reformulated 
with the variables Ha^ Ha^ and Tq, rather than E^ T, 
and kx. This mapping is beyond the scope of the present 
work, but we will outline how it can be done. 

The threshold applied magnetic field can be found by 
mapping kx = 7t/21x where lx{Ha) is the fiux penetra- 
tion depth from the Bean model. At the same time 
E = E{Ha,Ha).^^ For low E, we have T ^ Tq. The 
threshold temperature can be defined as the temperature 
where /^^ = 1 , i.e., full penetration is reached without 
the instability being nucleated. In this limit the domi- 
nant mechanism for prevention of the instability is the 
heat transfer to the substrate. Hence, a simple approxi- 
mation for the threshold temperature can be obtained by 
kx -^ 0. Then ah three cases, [Eqs. ([26|, ([32|, and ([35|], 
give the same condition: E = P/nF. 



B. The stability diagram 

Figure [3] shows a linear stability diagram in the T — E 

plane, calculated from Eq. (|23|), with parameters a — 

-5 



10-^ p = 0.1, 7 = 10, ni = 20, and Ix = 0.1. 
White color corresponds to values of T and E where 
max{ReA(k)} < 0, i.e., in the white regions the dis- 
tributions are stable. In the figure, there are two stable 
regimes. The lower one is the low-dissipative fiux creep 
state, which is stable for E < Etii{T). The upper is 
the high-dissipative fiux-fiow regime, which is stable for 
E > Jc{T)^ i.e., when the J—E curve is linear, see Eq. (|3|. 
These two stable domains are separated by an unstable 
region, where the arrow is meant to remind us that for 
the unstable state, E and T are bound to increase with 
time. 

The color in the diagram describes the proprieties of 
the most unstable mode: red is uniform {ky = 0), while 




FIG. 3. (color online) Stability properties of the supercon- 
ductor in the T — E plane. White denotes stable, coloured 
unstable. Red, green and yellow, means uniform oscillatory, 
non-uniform oscillatory, and fingering instability, respectively. 



The solid curves are E\ 



(0) 



r(i) 



and E, 



(2) 



Parameters are 



10" 



/3 = 0.1, 7=10, 



0.1, m = 20. 



green and yellow are non-uniform {ky > 0). Both the 
green and the red are modes giving oscillations (Im A ^ 
0). The colored lines correspond to the Eqs. (26), (32) 



and (35) derived, respectively, for limiting cases 0, 1 and 



2. As expected from the conditions of derivation, E^J (T) 

and E^ (T) follow nicely the edge of instability at low 
T. This means that at low T, the nucleated instability 
should be uniform, or close to uniform. At high enough 
T, the most unstable mode is at finite wavelength, i.e., it 
should giver rise to fingering structures. As expected, the 
instability threshold in this case is nicely approximated 

Inside the instability region, we have not derived any 
analytical expressions, since the linearized equations are 
not valid far from the instability threshold. However, it 
is possible to extract some qualitative information from 
the diagram also in this case. Of particular interest is 
that most of the diagram is yellow, which indicates that 
modes with ky > will grow fastest after the instability 
has been nucleated. 

The green line encloses the part of the diagram where 
the most unstable mode has ImA 7^ 0. i.e., it might be 
possible to detect damped oscillations. 



C. Simulations 

The simulations were carried out in two steps, as de- 
scribed in Sec. [llj First, the background state was pre- 
pared with thermal feedback turned off, giving uniform 
T = Tq. Second, the thermal feedback was turned on, 
and the perturbation 6T started evolving. 

Figure [4] shows successive T-maps at times ti < t2 < 
ts < ^4 after the thermal feedback was turned on, for 
two separate runs with Tq = 0.15 (left) and Tq = 0.4 
(right). We used the following values of the parameters: 
7 = 10-^ /3 = 0.1, 7 = 10, /^ = 0.2, and rii = 20. 
The ramp rates were chosen as Ha = 4 • 10~^ and 5 • 
10~^ for To = 0.15 and 0.4, respectively. Both values are 
just above the instability thresholds found heuristically 
by varying Ha- 

At ti = 10 the temperature is elevated in the flux- 
penetrated region. Even though the temperature rise is 
rather small, the instability is already nucleated and the 
appearance of a dendritic flux avalanche inevitable. For 
To = 0.15 the heated region is a narrow band near the 
edges, while for To = 0.4 it is much wider and extends 
almost to the flux front. 

The panels at ^2 show the T-maps when the fiux-flow 
hot spots appear. Each panel has only one hot spot, 
which is characterized by having the highest tempera- 
ture. At this early time it is still not much higher than 
the surroundings, and hence just barely visible. The 
hot spots appear because Jc decreases faster than J and 
eventually some position reaches the flux-fiow condition 
J = Jc The two runs develop at different rates, so we 
have t2 = 28 for To = 0.15 and ^2 = 12.5 for To = 0.4. 
Because the hot spots are characterized by the high flux- 
flow resistivity, they will quickly heat to the superconduc- 
tivity transition temperature. The locations of the hot 
spots are random, due to the uniformity of the sample. 

At ^3 = ^2 + 0.5 the avalanches have reached the prop- 
agation stage, where the hot spots have transformed to 
thin fingers. The fingers are either in the flux- flow or 
normal phase, and due to the high dissipation character- 
ized by extremely rapid propagation. The propagation 
is driven by the tip being adiabatically converted from 
the critical or Meissner state to the flux- flow state. At 
this stage, the propagation speed is not limited by ther- 
mal effects, so the speed of the front can even exceed the 
sound velocity. ^^ ^^ 

The final frames at ^4 = ^2 + 10 show the large branch- 
ing structures. The avalanches have reached their full 
extent and the structures are about to disappear as the 
heat is absorbed by the substrate. The dendritic struc- 
tures will remain in Bz and J^ 

Figure [4] allows us to present a fairly complete picture 
of how dendritic flux avalanches are nucleated and how 
they evolve. The avalanche has two distinct stages. In 
the first stage, the dynamics is characterized by the ther- 
momagnetic instability driven by a nonlinear I—V curve. 
Even though both E and T increase in time, they remain 
quite uniform. In the second - propagation - stage the 




FIG. 4. (color online) Evolution of dendritic flux avalanches, with To = 0.15 (left) and To = 0.4 (right). The temperature 
maps are at times ti < t2 < ts < t4, simulated with parameters a — 10~^, /3 = 0.1, 7 = 10, /a; = 0.2, rii — 20. 



dendritic flux structure is created. The dynamics is now 
totally driven by the highly dissipative (either flux-fiow 
or normal-state) branching structure invading the inner 
superconducting part of the sample. 

It is worth paying attention to the time spent in the 
two stages. The avalanche at Tq = 0.15 takes 28 time 
units in its first part to increase the temperature up to 
T = 0.16. In the propagation stage, it rises to T ~ 1 in 
less than 0.5 time units, and the creation of the branching 
structure takes 10 time units. This means that the actual 
instant of nucleation of the avalanche is much earlier than 
its first unambiguous signatures, such as the hot spot, or 
high pulses of T and E. 

Non-locality of the electrodynamics plays an essential 
role in the transition between the stages because it al- 
lows a uniformly nucleated instability to develop into a 
non-uniform avalanche. This behavior is different from 
parallel geometry, where the interaction between modes 
is absent, and a uniformly nucleated instability typically 
develops into a global flux jump.^ 



D. Comparison of the results 

Let us now compare the results of the simulations 
with the predictions of the linear stability analysis, us- 



ing identical parameters. Figure [5] presents a stability 
diagram in the T—E plane. Again, white regions are sta- 
ble and yellow ones are unstable according to the linear 
stability analysis. Here we do not distinguish different 
kinds of instability. The dots in the figure correspond 
to {Tmax, £^max} pairs extracted from the runs shown in 
Fig. |4] For each time, we collect the maximum temper- 
ature Tmax and the electric field T^max at the same time 
and position. Inside the instability region the points are 
ordered in time since Tmax(^+At) > Tmax(^)- This means 
that we can follow the evolution of the instability as a tra- 
jectory in the T— T^ diagram. The numbers correspond 
to the panels in Fig. |4] and link subsequent stages of the 
development to the temperature distributions. 

The onset of instability is the £^max with the lowest 
value. At Tq = 0.15 the lowest T^max is very close to 
the instability threshold calculated by the linear stability 
analysis. At Tq = 0.4 the lowest ^max is somewhat higher 
than the threshold. Hence we can say the linear stability 
analysis gives a good, but conservative, estimate for the 
actual instability onset. 

When the avalanche reaches the propagation stage, 
the trajectory {Tmax, £^max} makes a strong turn since 
^max ^ Jc{To) has reached it maximum value while T^ax 
increases with even faster rate. The maximum tempera- 
ture (not shown in the figure) is reached when the heat 



O 




FIG. 5. (color online) The development of an avalanche 
compared with the stability properties. Yellow means linearly 
unstable and white means stable. The points {Tmax,^max} 
are extracted from the numerical simulations. The blue points 
are for To = 0.15, red for To = 0.4. The labels 1-4 refer to the 
temperature distributions of Fig.[4| Parameters are a — 10~^, 
13 = 0.1, 7 = 10, U = 0.2, m = 20. 



ble region shrinks significantly when (Jra is increased from 
to 100. As expected, the metal layer mainly affects the 
stability at high E and T, when the conductivities of the 
two layers are within the same order of magnitude. 

Of particular importance is that Fig. [6] predicts a 
threshold temperature Ti, indicated by dotted vertical 
lines in the figure. Above this temperature the system 
is stable no matter the value of E. This opens a possi- 
bility that an avalanche can terminate without heating 
the sample above Tc. Thus the bistable properties of 
the system depend crucially on the value of Ti, which, 
in turn, depends on cFm- For example, Gm — 1 gives 
T\ = 0.97. In this case, a dendritic flux avalanche will 
most likely develop in the same way as without a metal 
layer. For am = 10, we have Ti = 0.8, which gives some 
prospects of avalanches terminating without the super- 
conductor being heated above Tc. For the highest con- 
ductivity, am = 100, we have Ti = 0.5, which means 
that it is likely that the magnetic braking will suppress 
formation of dendritic flux structures. 

Note that Ti is fundamentally different from the 
threshold temperature, Tth, often observed experimen- 
tally.^^ In particular, the threshold described by Tth is a 
consequence of the rapid growth of the thermal param- 
eters. It depends on E^ and does not alter the bistable 
properties of the system. 



removal is able to balance the heat production. After 
that, the temperature and electric field decrease with rel- 
atively slow rates until T^ax = 1 and T^max <C 1. The 
system is then again in the stable flux creep state and 
the avalanche is over. Then T^max will drop to a very 
low value, T^max ^ Hal^ diiid Tmax will dccrcasc with a 
rate determined by the Newton cooling in Eq. ([6| until 

J- max -^ D • 



V. METAL-COATED SUPERCONDUCTOR 

We will now consider how the thermomagnetic stability 
is affected by an adjacent normal-metal film. As depicted 
in Fig. [l] the metal layer is close, but not in a thermal 
contact with the superconductor. The linear stability 
analysis of this model was done in Sec. |TI| 



A. Dependence on the metal conductivity 

Let us first consider what happens when changing the 
normal metal conductivity am- Figure [6] shows the con- 
tours max{ReA(k)} = for a = 10-^ /3 = 0.1, 7 = 10, 
Ix = 0.1^ rii = 20. The contours are calculated by nu- 
merical solution of Eq. ([23]); the curves correspond to 
different am- The figure shows that the size of the unsta- 
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FIG. 6. (color online) The effect of changing the normal 
metal conductivity. The stability threshold contours in the 
T—E plane, for am — 0, 1, 10, and 100. Increasing metal layer 
conductivity improves stability at high E and T. Parameters 
are a = 10"^ /3 = 0.1, 7 = 10, l^ = 0.1, rii = 20. 
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FIG. 7. (color online) Stability properties of the metal- 
coated superconductor in the T — E plane. White denotes 
stable, coloured unstable. Red, green and yellow, means uni- 
form oscillatory, non-uniform oscillatory, and fingering insta- 
bility, respectively. The solid curve is E^'^^K Parameters are 
Gm = 100, a = 10-^ 13 = 0.1, 7 = 10, l^, = 0.1, m = 20. 



The offset of instability is thus appearing at high electric 
fields, of the order oi E r^ Jc/(^m- 



C. The stability diagram 

Figure [7| shows a T— £^ stability diagram calculated by 
numerical solution of Eq. (23) with cjm = 100. Other 
parameters are the same as in Fig. [3J As in the previous 
plot, white is stable, red is uniform oscillatory, green is 
nonuniform oscillatory, and yellow is fingering instability. 
For low T and E, the result is just as for the uncoated 
sample. Consequently, the low-T threshold conditions 



.(0) 



(1) 



for onset of instability, E^^i ^^^ ^th ' should be valid 
also for the metal-covered sample. For high T or E^ the 
differences compared to the uncoated sample are sub- 
stantial. According to the diagram, at £^ > Jc/cTm and 
T > Ti ~ 0.5 the system is stable. 

The upper edge of the instability region is offset of the 
instability, where all modes become stable. Yet, many 
modes will have values Re A ~ which means that they 
are almost stationary. We thus expect that avalanches 
subjected to magnetic braking will give rise to oscillations 
with long lifetimes in electric field and temperature. 

The analytical curve, Eq. ( |38| ), derived on the assump- 
tion of long wavelengths, provides a very good fit for the 
instability offset threshold, except for the temperatures 
close to Ti. 



B. Threshold for offset of the instability 



D. Simulations 



Let us now derive analytical expressions for the con- 
ditions for offset of the instability at high E and T. We 
assume that E and T are constant solutions of the non- 



linear equations, so that Eq. (22) can be used as a formal 



solution of the linearized equations in Fourier space. At 
the high electric fields, the background state may have 
evolved significantly from what is described by the Bean 
model, and the solution can be non-stationary, typically 
with T > To and £; > Ha. 

The magnetic braking comes into play when the non- 
linear exponent of the composite system, n from Eq. (15), 
is reduced. When nsJm ^ Js, we have 



1 + Js/Jn 



(36) 



where Jm = cr^y^E and Jg ^ Jc-, when Ug^ I. 

Several modes may be important for the offset of insta- 
bility. We will here focus on uniform oscillatory modes, 
which can be found by solving Re A = with ky = 0. 
Hence, 

B = kl^^ {akl + p) dm-Y (Jc - cJmE) F = 0. (37) 



Solving for E gives 






akl + /3 



(38) 



Figure [8] shows temperature distributions at times ti < 
t2 < ts < t^ < ^5 obtained by numerical simulations. 
Except the presence of the metal layer with dm = 100, 
the calculation procedure and the initial parameters are 
identical to those used in Fig. |4] The times ti to t^ are 
the same as in Fig. [4j while the last frame is at a much 
later time, ^5=^2 + 30. The time evolution of the two 
runs with Tq = 0.15 and Tq = 0.4 are quite different, so 
we discuss them separately. 

For To = 0.15, the states at ti and ^2 are exactly the 
same as for the previous run with the uncoated sample. 
This means initial phase of the avalanche is not affected 
by the metal layer. At ^3, a hot spot heated to T ^ 0.25Tc 
is visible. It is propagating away from the edge, but much 
more slowly than the propagating finger of Fig. [4j At ^4, 
the state has deviated further from that for the uncoated 
sample and there are even two fingers developing in par- 
allel. The temperature inside the fingers is approximately 
T ~ 0.5. The final frame shows the temperature distri- 
bution at ^5 = ^2 + 30, which is more or less the final 
development of the instability. The two structures have 
at this point exceeded the flux front but they are still 
much smaller than the dendritic flux avalanche in Fig. [4] 

For Tq = 0.4, the temperature distribution at nucle- 
ation is the same as in the uncoated sample. However, 
all times t > ti are affected by the magnetic braking. 
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FIG. 8. (color online) Evolution of flux avalanches subjected to magnetic braking, with To = 0.15 (left) and To = 0.4 (right). 



The temperature maps are at times ti < t2 < ts < t^ < t^, simulated with parameters a = 10 , /3 = 0.1, 7 = 10, 



0.2, 



ni 



20. 



In the figure, it is impossible to identify one individual 
avalanche. Instead, there is a quasi-periodic structure, in 
which the peaks are heated to T ~ 0.56. The structures 
are mainly in the critical state region and do not pen- 
etrate into the Meissner state region. In the last frame 
taken at ^5=^2 + 30 the temperature is almost back to 
To- 

The comparison between Figs. |4] and [S] show that, for 
To < Ti and E <C Jd^m^ the initial stage of the insta- 
bility development is not affected by the presence of a 
metal layer. However, with high metal layer conductiv- 
ity, the magnetic braking prevents the rapid propagation 
of avalanches, and as a result the final structure of the 
avalanche is strongly altered. The avalanche at Tq = 0.15 
looks like a blob, with some tendency of branching, while 
at To = 0.4 it becomes a semi-periodic wiggling of the 
fiux front. In both cases, the magnetic braking prevents 
formation of the dendritic structure. 



E. Comparison of the results 

In order to check the accuracy of the predictions of the 
linear stability analysis, let us compare them with the 
results of simulation using identical parameters. Figure |9] 
shows di T — E stability diagram, where the background 



color represents results of the linear stability analysis, so 
that white means stability, yellow instability. The points 
are {Tmax, £^max} pairs extracted from the simulations, 
^max is the highest temperature at a given time and T^max 
is the electric field at the same time and location. The 
blue points come from the run nucleated at Tq = 0.15, 
red points at Tq = 0.4. The numbers 1-5 correspond to 
the panels in Fig. [SJ 

At To = 0.15, the instability is nucleated at low E and 
between ti to ^2, both T and E grow with time, exactly 
as without metal coating. At ^3, the electric field satu- 
rates at T^ ~ Jc/o'm^ Eq. (38), which is much lower than 
the saturation level T^ ~ 1 in the pristine superconduc- 
tor of Fig. [5J The fit with the linear stability analysis is 
remarkable as {Tmax,^max} follows the edge of the dia- 
gram. At t > ^5, the electric field and temperature drop, 
and the avalanche is over. This means that the simu- 
lated avalanche to large extent behaves as predicted by 
the linear stability analysis. 

The temperature Tq = 0.4 is close to Ti ~ 0.55. 
Hence, the {Tmax^^max} points are close to the insta- 
bility threshold at all times, and the trajectory devi- 
ates strongly from that of the pristine superconductor 
of Fig. [5] The electric field grows until it hits the max- 
imum at T^ ~ Jcl^rai whereupon the electric field stays 
constant for some time while T increases. At ^4, the 
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FIG. 9. (color online) The development of avalanches 

compared with the stability properties for a metal-coated 
sample. Yellow means unstable, white stable. The points 
{T'maxj^max} are cxtractcd from simulations: the blue points 
are for To = 0.15, red for To = 0.4. The labels 1-5 refer 
to the temperature distributions of Fig. [S] Parameters are 
(Tm = 100, a = 10"^ P = 0.1, 7 = 10, h = 0.2, m = 20. 



avalanche has reached its largest extent and the temper- 
ature and electric field decrease rapidly. At time t^ the 
electric field is so low that it is out of the range of the 
figure. 

Because the avalanches at the propagation stages fol- 
low the edges of the instability region, we can conclude 
that the linearized theory describes properly also the off- 
set of instability due to the magnetic braking. 



VI. DISCUSSION AND SUMMARY 

Type II superconducting films in the critical state 
are susceptible for the thermomagnetic instability, which 
may cause dendritic flux avalanches. Such avalanches 
are potentially damaging for applications. Consequently, 
it is desirable to work out criteria telling when the in- 
stability appears and how its impact can be minimized. 
Derivation and analysis of those criteria for supercon- 
ducting films (both uncoated and coated by a normal 
metal) was the purpose of the present work. To reach 
this goal we have performed both linear stability analy- 
sis and numerical simulations of the equations governing 
onset and propagation of coupled fluctuations of mag- 
netic fiux and temperature. Comparison of these results 
reveal the physical picture of the flux avalanche dynamics 
in thin-film superconductors. 



The state prior to linearization was quantified by the 
electric field E and temperature T. Perturbations of this 
state were analyzed in the Fourier space, their scales hav- 
ing been quantified by wave vector k. We considered the 
most unstable modes with k corresponding to the largest 
instability increment, max{ReA(k)}. In x direction, the 
most unstable mode is always the one with shortest wave- 
length, and thus the most unstable mode was associated 
with the flux penetration depth l^ = 7r/2kx. The ma- 
terial properties were combined into the dimensionless 
heat conductivity a, the dimensionless coefficient of heat 
transfer to the substrate /3, the Joule heating parame- 
ter 7, and the flux creep exponent parameter rii. These 
quantities depend only on measurable parameters, there- 
fore the theory is quantitative. The derived expressions 
were explored under assumption of realistic, for conven- 
tional superconductors and MgB2, temperature depen- 
dencies of the material parameters. 

The linear stability analysis showed that the super- 
conductor was stable for small E and Ix^ and high T, 
in agreement with previous results. Increased 7 and ni 
decreased stability and increased a and /3 improved sta- 
bility. Analytical expression where derived in three cases: 
(i) uniform oscillatory instability for small T, E^ and Ix] 
(ii) nonuniform oscillatory instability for intermediate T 
and E and for large Ix] (iii) fingering instability for higher 
T and E. 

We have solved the same set of equations numerically, 
but without linearization and with more realistic bound- 
ary conditions. The initial state was numerically pre- 
pared by ramping of the applied magnetic field, with 
thermal feedback tuned off. At a given flux penetra- 
tion depth, the thermal feedback was turned on, and the 
evolution of the perturbation was explored. We report 
the results for parameter combinations just above the 
heuristically found instability threshold, where the insta- 
bility developed into a dendritic flux avalanche. We have 
found that the avalanche has two distinct stages. First, 
just after the nucleation, the sample is in the flux creep 
state, and even though T and E increased exponentially 
with time, the perturbations are almost uniform. Since 
E is low, the flux traffic is limited. Second, the prop- 
agation stage begins with the appearance of a hot spot 
in the flux- flow phase. The hot spot soon turns into a 
thin finger, which rapidly propagates into the sample, 
creating a branching structure. Hence, we conclude that 
the avalanche is initiated by the thermomagnetic instabil- 
ity, but this is not the main mechanism behind creation 
of the dendritic structure. Instead, it is more accurate 
to think about the dendritic structure as created by a 
highly-dissipative phase (flux-flow or normal) invading a 
low-dissipative one (flux creep or Meissner). 

For two different temperatures, the instability onset 
found heuristically based on numerical simulations was at 
slightly higher electric fields than those predicted by the 
linearized theory. Therefore, the linear stability analysis 
should provide a good, but conservative, estimate of the 
instability threshold. 
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The magnetic braking was described by a simple model 
where the superconducting and metallic layers were con- 
sidered as two electrically and thermally isolated current- 
carrying layers connected in parallel. 

Linearization of the set of equations for the coated sys- 
tem gave two new parameters: the normal metal sheet 
conductivity am and the nonlinearity exponent of the 
composite system n = n{E^T). For low T and E the 
conditions for onset of instability where unchanged by 
the presence of the metal layer. At higher E, when 
TisJm ^ Jc^ the limit for offset of instability was lowered 
by the presence of the metal. The effect was stronger 
for increasing a^. The theory predicts that there is a 
temperature Ti{am) < Tc, such that the system is stable 
for any E when T > Ti. This means that the system 
can recover from an instability without being heated to 
the normal state. Analytical expressions for offset of the 
uniform oscillatory instability by magnetic braking were 
derived in the low temperature limit. 

Numerical simulations confirmed the prediction of the 
linear stability analysis that the impact of avalanches 
may be significantly reduced by magnetic braking. The 
avalanches subjected to magnetic braking created ex- 
tended protrusions of the flux front rather than dendritic 
structures. The reason for the suppression of the flux 
traffic is that finite conductivity of the coating layer lim- 
its the electric field. The maximum temperature during 
the avalanches in a coated film did not exceed Tc^ just as 
predicted by the linearized theory. 

The main conclusions from this work are that a uni- 
form thermomagnetic instability can develop into a den- 
dritic flux avalanche. In the long, initial phase of the 
avalanche, the electromagnetic non- locality causes the 
appearance of a small non-uniformity in the tempera- 
ture perturbation. A dendritic structure is created when 
a hot spot appears in a random position, and develops 
into a finger, which propagates away from the edge, at 
very high velocity. The magnetic braking does not affect 
the nucleation of the thermomagnetic instability at low 



electric fields, but it may significantly reduce the impact 
of avalanches, and may suppress formation of dendritic 
structures. 

Several predictions from this work can be checked ex- 
perimentally. The prediction that the instability at low 
temperatures is nucleated uniformly along the edge is 
not easy to check directly, since it is very difficult to dis- 
tinguish a stable configuration from an unstable one in 
avalanches at early times. Instead, one can look for col- 
lective oscillations, e.g., by using an array of Hall probes 
along the edge, since oscillations are the hallmarks of 
modes with ky < kx- The formulas for the instability 
threshold could be checked indirectly, by finding and fit- 
ting values for the threshold length lx{Ha) with exper- 
imental values. To check the theory for offset of insta- 
bility by magnetic braking one could also use an array 
of Hall probes, this time at the edge of a superconduct- 
ing strip covered by metal. The electric field observed 
during the avalanche should be inversely proportional to 
the sheet conductivity of the normal metal. The frequen- 
cies uj = yJCjA of the undamped oscillations also bear 
information of physical parameters of the system. 

The theory can be modified to handle the case when 
the metal layer and superconductor are in a close thermal 
contact. In this limit, the thermal shunt should shift the 
onset of instability, because normal metals have k^ oc T 
and Cm ex T at low temperatures. The magnetic braking 
should, however, be less efficient than the case considered 
in this work, since the Joule heating in the normal metal 
will also heat the superconductor. Therefore, the Joule 
heating term should be enlarged from 2s^ to j£^. In 
the general case on must model the temperature in the 
superconductor and normal metal separately. 
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